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ABSTRACT 

To investigate the origins of SO galaxies, we present a new method of analyzing their 
stellar kinematics from discrete tracers such as planetary nebulae. This method in- 
volves binning the data in the radial direction so as to extract the most general possible 
non-parametric kinematic profiles, and using a maximum likelihood fit within each bin 
in order to make full use of the information in the discrete kinematic tracers. Both 
disk and spheroid kinematic components are fitted, with a two-dimensional decompo- 
sition of imaging data used to attribute to each tracer a probability of membership in 
the separate components. Likelihood clipping also allows us to identify objects whose 
properties are not consistent with the adopted model, rendering the technique robust 
against contaminants and able to identify additional kinematic features. 

The method is first tested on an N-body simulated galaxy to assess possible sources 
of systematic error associated with the structural and kinematic decomposition, which 
are found to be small. It is then applied to the SO system NGC 102 3, for which a plan- 
etary nebula catalogue has already been released and analyzed by iNoordermeer et al.l 
(|2008l ). The correct inclusion of the spheroidal component allows us to show that, 
contrary to previous claims, the stellar kinematics of this galaxy are indistinguishable 
from those of a normal spiral galaxy, indicating that it may have evolved directly from 
such a system via gas stripping or secular evolution. The method also successfully 
identifies a population of outliers whose kinematics are different from those of the 
main galaxy; these objects can be identified with a stellar stream associated with the 
companion galaxy NGC 1023A. 

Key words: galaxies: elliptical and lenticular - galaxies: evolution - galaxies: in- 
dividual: NGC 1023 - galaxies: individual: NGC 1023A - galaxies: kinematics and 
dynamics. 



1 INTRODUCTION 

Lenticular or SO galaxies lie between ellipticals and spirals 
on the Hubble sequence, since they have the featureless old 
stellar populations of elliptical systems, but also contain the 
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disk components associated with spirals. They thus poten- 
tially represent a key element in attempts to understand 
the relationship between the main types of galaxy, but at 
the moment there is no consensus as to the end of the Hub- 
ble sequence to which they are most closely related. Clearly, 
some process has consumed their gas and shut off their star 
formation. However, this termination could be the result of 
galaxy mergers, much in the manner that star formation is 
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believed to be quenched in elliptical galaxies, or it could 
arise from some much gentler process, such as ram pressure 
stripping, which simply removes the gas from normal spiral 
galaxies. 

Perhaps the best way of discriminating between these 
scenarios is offered by the stellar dynamics of SO galaxy 
disks. If they are the product of relatively gentle gas- 
stripping processes, one would expect the stellar dynam- 
ics to be unaffected, and should be identical to those of 
the progenitor spiral system, with rotation dominating rel- 
atively modest amounts of random motion in the disk 
ljArag6n-SaIamancall2008l : IWilliams et al.ll2010h . If, however, 
the SOs are the product of minor mergers, with a mass ratio 
higher than 10:1, one would expect the merging process to 
be imprinted in the stellar dynamics, with random velocities 
as hi gh or even higher than the rotational velocities at all 
radii (iBournaud et al.ll2005l '). 

One difficulty in using such a diagnostic is that it re- 
quires accurate stellar kinematics measured in the outer 
parts of the galaxy where the disk dominates, and at 
low surface brightness, making conventional absorption-line 
spectroscopy difficult. However, planetary nebulae (PNe) 
provide an excellent probe of stellar kinematics in this 
regime: they are easily detected and can have their ve- 
locities measured from their characteristic strong emission 
lines, and they have proved t o be an unbiased tracer of 
the global stellar population (ICiardullo et al.1 1 19891 . |2002| : 
iNapolitano et al.l l200ll : ICoccato et all 120091 '). We have de- 
signed and built a special-pur pose instrument, the Planetary 
Nebula Spectrograph [PN.S: lOouglas et al.1 l|2002l )]. specifi- 
cally to exploit this tracer of stellar kinematics at low surface 
brightnesses . Although originally intended to study ellipti- 
cal galaxies jRomanowskv et al.1 l2003l : ICoccato et al.l l2009l : 
INapolitano et al.l 20091 '). PN.S has already pro ved very effec- 
tive i n exploring the disks of spiral galaxies (|Merrett et al.l 
I2006D . so an obvious next step was to use it to probe the 
kinematics of SO systems. 

As a pilot study, we observed the relatively nearby SO 
galaxy NGC 1023 with PN.S. The observations and the re- 
sulting catalogue of positions and line-q f -sight velocities of 
PNe are described in lNoordermeer et al.l l|2008l ). The analy- 
sis in that paper concluded that NGC 1023 has very peculiar 
kinematics, giving the appearance of a normal rotationally- 
supported disk galaxy at small radii, but having entirely ran- 
dom velocities at large radii, inconsistent with either of the 
expected scenarios. However, the relatively simple dynami- 
cal analysis that led to these conclusions had some signifi- 
cant shortcomings. First, it assumed that the light at large 
radii was dominated by the disk component, and neglected 
any contribution from the spheroid. As we will see, this as- 
sumption can lead to sizable systematic errors in the inferred 
disk kinematics. Second, it calculated dynamical quantities 
such as velocity dispersions by binning the data both az- 
imuthally and radially in the galaxy. Although such binning 
provides an adequate signal from a sparsely-sampled veloc- 
ity field and a non-parametric measure of the local kine- 
matics, as we will show below the kinematic properties vary 
continuously and quite rapidly with azimuth, so such bin- 
ning averages away significant amount of information that is 
present in the raw data. The binning process also makes it 
difficult to identify any contamination from PNe that are ei- 



ther unrelated background objects, or lie within the system 
but do not match the expected kinematics. 

In this paper, we therefore present a somewhat more so- 
phisticated analysis designed to circumvent these shortcom- 
ings. While still binning the data in the radial direction, so as 
to extract a general non-parametric view of the galaxy's dy- 
namics, we apply a maximum likelihood analysis within each 
radial bin, so as to extract the maximum amount of informa- 
tion from the azimuthal variations in kinematic properties. 
We also model both disk and spheroid in the system, and 
use a full two-dimensional photometric decomposition of the 
galaxy, which allows us to allocate to each PN a probability 
as to the component to which it belongs. 

The method we have developed to construct such a 
model is described in Section [21 and in Section [3] we test 
the technique on a simulated galaxy for which we have 
prior knowledge of the structural components and kinemat- 
ics. In Section U we make a first real application of the 
method to the NGC 1023 catalogue, in which we show that 
this likelihood fitti n g can reproduce the strange results in 
iNoordermeer et al.1 (|2008l ) if we also neglect the spheroidal 
component, but that the more complete model presented 
here results in a picture in which NGC 1023 has a normal 
rotationally-supported disk. We also illustrate the power of 
likelihood fitting by showing how outlier points can be iden- 
tified in the catalogue, and associated with a minor on-going 
merger. The conclusions are presented in Section [S] 



2 KINEMATIC LIKELIHOOD FITTING 

In this section we give a brief introduction to likelihood 
analysis, before going on to the specifics of fitting partic- 
ular galaxy models. In Section 12.11 we present the case of 
a pure disk galaxy, as it provides both a simple example of 
the technique and a useful point of comparison to previous 
analyses that have assumed that the disk is the dominant 
component. In Section [2.21 we proceed to the more realistic 
situation in which we model both disk and spheroidal com- 
ponent. In Section r2.3l we discuss how the individual discrete 
kinematic tracers can be assigned, at least in a probabilistic 
sense, to either the disk or the spheroid using photometric 
data. 

In general, given a set of A*' independent velocity mea- 
surements, Vi = (ui, ujv), drawn from a probability den- 
sity function F{vi; 0), where 6 = {V, a) is a set of parameters 
whose value is unknown, the likelihood function is given by 



c = \{F{vi-e). 



(1) 



The values of V and a that maximize C are the best estima- 
tors for the true values of these parameters. Moreover the 
method of maximum likelihood coincides with the method 
of the least squares in the special case of a set of N Gaus- 
sian distributed independent measurements, in which case 
the likelihood function is directly related to the usual 
statistic. 



Ax (f) = -2Aln£(6i). 



(2) 



Thus, in this case it is straightforward to determine the con- 
fidence region around the best estimator: 



ln£(6») ^ ln£™, - Aln£, 



(3) 
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where values of Ax^ or — 2Aln£, corresponding to a desired 
confidence limit, for joint estimation of m parameters, are 
readily available in tabulated form. In this paper we use 
a Ax^ that corresponds to a la coverage probability. In 
particular, Ax^ ~ 5.39 for m = 4, and Ax^ = 4.11, for 
m = 3. 

Finally, once the best-fit parameters have been found, 
one can calculate the contribution of each data point to the 
likelihood given this optimum model. To quantify this con- 
tribution in terms of whether the kinematic tracer in ques- 
tion is consistent with being drawn from the best-fit model, 
one can generate a large number of individual velocities 
drawn from a velocity distribution matching the model, and 
see at what percentile of the distribution the data point lies. 
We can thus identify any objects whose velocities lie outside 
this confidence interval, which are most likely interlopers. 
By rejecting these data points and iterating on the fit, we 
can render this process robust against a small amount of 
contamination, in a straightforward gene ralization of the 3(7 
clipping procedure (jMerrett et al.ll2003l l used in determin- 
ing membership of velocity distributions that are assumed 
to be Gaussian. We define the likelihood clipping probability 
threshold as the value at which we cut the distribution, and 
we reject all the PNe beyond this limit. 

In principle, such likelihood fitting could be applied to 
a full dynamical model across an entire galaxy, at all radii 
and azimuths. However, in this analysis we are interested in 
leaving as much freedom as possible in the resulting radial 
profile of dynamical properties, so as to explore the full range 
of possibilities predicted by different models of SO formation. 
We therefore adopt a hybrid approach in which the data are 
binned into elliptical annuli matched to the inclination of the 
disk, to extract points over a range of intrinsic radii. This 
binning allows each section of the disk to be modeled as 
an independent non-parametric data point. However in the 
azimuthal direction around each bin we use a full likelihood 
analysis that accurately models the expected variation in the 
line-of-sight velocity distribution with azimuth, as derived 
from very general dynamical considerations. 

2.1 Likelihood analysis for a disk model 

Consider a general rotating disk model for a galaxy. Part 
of the line-of-sight velocity of each object within it, Vi, is 
the projection of the galaxy's mean rotation velocity V , at 
an azimuthal angle cj> within the galaxy, corrected for the 
inclination of the galaxy to the line of sight i, and for the 
systemic velocity of the galaxy itself Vsys'- 



Vios = Vsya + V sin{i)cos{<f)). 



(4) 



Superimposed on the net rotational velocity are the ran- 
dom motions of the individual stars, which can be quanti- 
fied by their velocity dispersion in different directions. For an 
axisymmetric disk, these components are most naturally ex- 
pressed in cylindrical polar coordinates aligned with the axis 
of symmetry, (i?, <j), z), and the line-of-sight velocity disper- 
sion, aios, is made up from a projection of these components 
such that 



2 2"2-'2 2"2- 2 2 2- 

^loa = sin i sin (j) + sin i cos (/> + cos i 



(5) 



galaxy, like the ones we are focusing on, it does not project 
much into the line-of-sight velocity dispersion, so the value 
of (Tjos is dominated by the other two components. As can be 
seen from the above equation, crfo^ varies sinusoidally such 
that its value is set by (tr on the minor axis and on 
the major axis. Thus, by fitting the variation in aioa with 
azimuth one can determine the values of both of these com- 
ponents of the galaxy's intrinsic velocity dispersion. 

If we adopt the simplest possible model in which the 
line-of-sight velocity distribution is Gaussian at every point 
with a mean velocity of Vios and a dispersion of aios , we now 
have the requisite form for the probability density function. 



F{vi-V,ar,a^) oc exp 



[n—Vios{V)? 



(6) 



Generally speaking, in disk galaxies, is the smallest com- 
ponent of the velocity dispersion, and in a nearly edge-on 



The values V, 0-^,0-^ that maximize C are the best estimators 
for the true values of these kinematics parameters within 
each bin. 



2.2 Likelihood analysis for a disk + spheroid 
model 

In addition to the disk component, most systems also usually 
contain a spheroidal component comprising either a central 
bulge or an extended halo (or both). Indeed, one of the tra- 
ditional defining features of SO galaxies is that they often 
have rather prominent bulges. A kinematic model without 
such a component may therefore be significantly incomplete; 
worse, the amount of bulge light varies with azimuth around 
the galaxy, since close to the minor axis in a highly-inclined 
system the bulge is usually the dominant component, so one 
might expect the variation in velocity dispersion with az- 
imuth, used above to extract the different components of 
disk velocity dispersion, to be systematically distorted. 

One slight complication in adding in such a spheroidal 
component is that it does not have the same symmetry prop- 
erties as the disk component. Thus, the elliptical annuli that 
contain data from a small range in radii in the disk does not 
correspond to the same range in radii in the more spherical 
spheroidal component, but samples a larger range of radii 
both due to the difference in shape and the effects of in- 
tegration along the line of sight. However, the variation in 
kinematics with radius in such a hot component is generally 
rather slow, so the greater averaging in radius imposed by 
the choice of elliptical annuli should not have a major im- 
pact. Further, the validity of this assumption of slow vari- 
ation with radius can be tested a posteriori by seeing how 
the inferred kinematics of the spheroid change from bin to 
bin. 

Accordingly, we adopt the simplest possible model for 
the kinematics of the spheroidal component, in which the 
line-of-sight velocity distribution is assumed to be a Gaus- 
sian with zero mean velocity (relative to the galaxy). We 
have thus assumed that any rotation in the spheroidal com- 
ponent is negligible, although, as we shall see below, this 
assumption can also be relaxed. The velocity dispersion, as, 
is left as a free parameter to be modeled by the data. The 
only other parameter that we need to specify is the proba- 
bility that each individual tracer object at its observed lo- 
cation belongs to the spheroidal component, fi, so that the 
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full probability density function can be written 



3 APPLICATION TO A SIMULATED GALAXY 



fi 

F{vi;V,ar,a^,as) oc — exp 



+ ^exp 



{V^ - Vlosf 



2af 

LOS 



(7) 



where the velocity dispersions in the denominators of the 
term in front of each Gaussian ensure that this function 
is properly normalized. From this velocity distribution one 
can then construct a likelihood function for a particular 
set of model parameters using Equation [l] The values of 
V,<jR,<j^,Os that maximize C then provide the estimators 
for these parameters in the best-fit model. 



As a first test of this method, we perform the fit to a 
model galaxy obtained from a self-consistent N-body sim- 
ulation. The simulation was provided by Kanak Saha (pri- 
vate communication). This simulated galaxy is constructed 
using a nearly self-consistent bulge-disk-dark halo model 
(jKuiiken fc Dubinski|[l995l : IWidrow et "al]|2008l ') to mimic a 
typical lenticular galaxy. The value of Toomre Q is high 
enough to prevent strong 2-armed spirals from forming in 
the disk. The spheroid follows a Sersic profile with an index 
of 3.5, while the disk is exponential with a vertical sech^ 
profile. The dominant dark matter halo specifies the sys- 
tem's rotation curve. We chose to "observe" this galaxy at 
an inclination of 60 degrees. For the following tests, we set 
the likelihood clipping probability threshold to 2.1a". 



2.3 Spheroid— disk decomposition 

The only parameters that we have not yet specified are the 
fi values, which define the probability that a given object 
is in the spheroid. We cannot solve for these quantities di- 
rectly in the likelihood maximization since, as noted above, 
the relative contributions of spheroid and disk vary with az- 
imuth within a single bin, so each kinematic data point has 
its own individual value. Fortunately, we have additional 
information that has not yet been used. In particular, we 
can apply a two-dimension al galaxy fitting routine such as 
GALFIT l|Peng et al.ll2002l ) to imaging data in order to de- 
compose the starlight into spheroid and disk components. 
The fit to the light distribution then provides a direct esti- 
mate of the fraction of the starlight from each component 
at the location of any stellar tracer. 

This approach for determining the decomposition into 
bulge and disk components has the great advantage that the 
broad-band imaging data offer a much less sparse sampling 
of the stellar spatial distribution than the PNe, providing 
an intrinsically more accurate answer. Furthermore, there 
are significant selection effects in detecting PNe, as they 
are harder to find against the bright background of the in- 
ner parts of a galaxy, so their spatial arrangement is not 
an unbiased representation of the stellar distribution. Note, 
however, that this bias is a purely spatial one, in that all 
line-of-sight velocities are equally detectable, so their use in 
the kinematic analysis is not in any way compromised. 

The only slight subtlety in applying such analysis to 
PNe is that their number per unit stellar luminosity has 
been shown to vary systematica lly with the colour of the 
population (I CiarduUo et al.lll99ll ). with a lo wer PN density 
per u nit galaxy luminosity for redder objects (jBuzzoni et al.l 
I2OO6I ). If there is a difference in colour between the disk and 
spheroidal components, then one has to apply a correction 
in order to convert the fraction of spheroidal light at any 
given point into the probability that a PN detected at that 
point belongs to the spheroid. In practice, such colour terms 
can be straightforwardly determined by performing the de- 
composition on images taken in different bands, and using 
the resulting colours of the different components to correct 
the probability. 



3.1 Testing the likelihood method using a priori 
knowledge of PN positions 

In a simulated galaxy, we have access to inside informa- 
tion which allows us to test specific aspects of this mod- 
eling process. In particular, the "stars" are tagged accord- 
ing to whether they are members of the disk or spheroid, 
so we can assign the objects to specific components with 
certainty, rather than using the probabilistic approach de- 
scribed above. Accordingly, we have carried out the likeli- 
hood analysis without the probabilistic decomposition, with 
the results presented in Figure [T] The upper panel shows 
the results obtained with a generous catalogue of 437 ob- 
jects, while the lower panel shows how well we can do with 
only 136 kinematic tracers. Although not strictly valid for 
kinematically-hot SO galaxies, an asymmetric drift correc- 
tion has been applied to convert the derived rotation veloc- 
ity into a circular velocity for direct comparison with the 
simulation's known functional form, using the formula 

=V + a^ - a^ 1 h -71 , (8) 

\ td amr J 

where Vc is the circular velocity, is th e median radius i n 
each bin, and tb is the disk scale length (jKormendvlll984l l. 
The final gradient term of the equation has been estimated 
assuming that a^ follows an exponential profile and perform- 
ing a linear fit between Ina^ and the radius. 

While the errors are, as expected, larger for the smaller 
sample, there are no systematic differences between them, 
and both do a good job of recovering the simulation's cir- 
cular velocity. It is interesting to see that ajj is moderately 
but systematically too low for the small sample-size, and 
remains low at large radii even for the larger sample. Since 
the main remaining assumption in this model is that the 
velocity distributions of the individual components are in- 
trinsically Gaussian, it seems likely that this modest system- 
atic error occurs due to the breakdown of this assumption. 
However, we note that a 4, and v are both quite accurately 
recovered, so we can estimate with some confidence the bal- 
ance between random and rotational motion, o^/v, which is 
the main physical quantity that we are seeking to use as a 
diagnostic in this analysis of the origins of SO galaxies. 
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Figure 1. Derived circular velocity and tlie components of tiie 
velocity dispersion versus radius for the model galaxy in the disk 
and in the bulge, for 437 particles, upper panel, and 136 par- 
ticles, bottom panel. The filled symbols are from the maximum- 
likelihood analysis, with vertical error bars indicating uncertainty 
and horizontal error bars showing the extent of each radial bin 
(with the point plotted at the median radius for a PN in that 
bin). The dashed lines represent the model circular velocity, ve- 
locity dispersion in the disk and in the spheroid. 



3.2 Testing the effect of the photometric 
spheroid— disk decomposition 

We now introduce the remaining aspect of the full model- 
fitting process, the spheroid-disk decomposition, which al- 
lows us to assign a probabilistic membership of each kine- 
matic tracer to each component, since, as discussed above, 
in a real galaxy we do not have the luxury of this knowl- 
edge a prion. We constructed an image of the simulated 
galaxy and convolved it with a suitable Gaussian point- 
spread function. We then used GALFIT to model the result- 
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Figure 2. Results from the isophotal analysis of the simulated 
galaxy, black filled triangles, and the galaxy models obtained with 
GALFIT leaving the Sersic index n as a free parameter, red filled 
pentagons, imposing Sqalfit = 3.5, green open pentagons and 
imposing Sqalfit = 4, blue starred symbols. The top row shows 
the fitted position angles, the middle row the ellipticity and the 
bottom row the surface brightness. 



ing simulated broad-band image. The principal parameters 
returned by this fitting process are the spheroid and disk 
magnitudes, their respective scale lengths, the Sersic index 
of the spheroid, and the flattening and position angle of the 
components. 

Interestingly, if the parameters are all left free, then 
the spheroid is found to have a best-fit Sersic index of 2.8, 
significantly smaller than the known value for this simulation 
of 3.5. The value of the effective radius is also found to be 
systematically too small. However, an almost equally-good 
fit is found if we fix the Sersic index to 3.5. In fact, this 
fit is by some measures superior: Figure [2] shows the values 
of position angle, ellipticity and surface brightness obtained 
by fitting elliptical isophotes to both the simulated galaxy 
image and the models in which the Sersic index is either 
fixed or left free. While both models reproduce the position 
angle and surface brightness equally well, the ellipticity is 
clearly better fitted by the model with the Sersic index fixed 
at the right value. This conflicting information underlines 
the complexity of non-linear model fitting, and illustrates 
its basic limitations. 

Of course, for a real galaxy, we would not know the 
"right" value for the Sersic index, so would not be able 
to choose between these models. Such systematic errors in 
GALFIT fitting therefore pose a potential limitation to the 
effectiveness of the modeling procedure set out in this paper, 
if the resulting kinematics turn out to depend sensitively on 
the spheroid-disk decomposition. To assess the impact of 
such effects, we carried out the maximum likelihood mod- 
eling using both of these decompositions. We also checked 
the effect of kinematic tracer sample size by simulating cat- 
alogues of 437 and 136 PNe. The resulting kinematics are 
presented in Figure |3l The good news is that the results 
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Figure 4. Positions and velocities of the PNe detected in 
NGC 1023. Colors bar indicates velocities. The elliptical annuli 
show the radial binning of the data. The circled points are those 
rejected in the pure disk model, with the likelihood clipping prob- 
ability threshold increased to 1.65cr. Square s show PNe associ- 
ated with the companion galaxy NGC 1023A l|Noordermeer et al.l 
l2008f) . 



are largely insensitive to systematic errors arising from the 
disk - spheroid decomposition. For the smaller catalogue, 
the extra uncertainty that arises from the decomposition 
process mean that the errors in the derived kinematic quan- 
tities become quite large, but there is no evidence for any 
systematic error in the process. Thus, it appears that this 
maximum likelihood procedure is quite reliable, and robust 
against the most likely sources of systematic error in the 
spheroid-bulge decomposition. 



4 APPLICATION TO NGC 1023 

Having developed this methodology for extracting kinematic 
properties of disk galaxies from individual stellar tracers, 
and tested its reliability, we can now apply it to the ex- 
isting PNe data for the SO galaxy NGC 1023. The cat- 
alogue of 183 PNe posit ions and velocities, published by 
iNoordermeer et all (|2008l ). are presented in Figure |4] Note 
the "hole" in the middle of the galaxy resulting from the 
difficulty of detecting PNe against the bright stellar conti- 
nuum at these radii, as discussed above. Application of the 
new method to these data has the benefit that they have 
alread y been studied using a rn ore conventional tilted-ring 
model l|Noordermeer et ahlbOOSl l , with which our results can 
be compared, plus the somewhat peculiar kinematics ap- 
parently found in this system warrant further investigat ion. 
Since the previous analysis by INoordermeer et al.l (|2008l ) ne- 
glected any contribution to the kinematics from the bulge, 
we begin by considering the disk-only model to make a direct 
comparison. 
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Figure 5. Derived mean rotation speed and the components of 
the velocity dispersion versus radius for NGC 1023, for the case 
of a disk-only model. The filled symbols are from the maximum- 
likelihood analysis, with vertical error bars indicating uncertainty 
and horizontal error bars showing the extent of each radial bin 
(with the point plotted at the median location for a PN in that 
bin). The open symbols reproduce t he results of the tilted-ring 
analysis by INoordermeer et al.l l|200^ 'l . 



4.1 The disk model 

The application of the above approach, in which the data 
are binned radially (into the annuli indicated in Figure |4| 
but modeled azimuthally by likelihood analysis, is presented 
in Figure [5] For this analysis and the following sections, we 
have estimated the galaxy's inclination from the ellipticity 
of the disk component derived from the two-dimensional fit 
to the photometry (see Section [4. 2p : assuming that the disk 
is intrinsically axisymmetric, we infer a value of i = 74.3 
degrees. Likelihood clipping has been applied such that 
PNe with a probability of less than 2.\a of being drawn 
from the disk model are rejected. As is clear from this fig- 
ure, the analysis largely reproduces the peculiar result of 
INoordermeer et al.l l|2008l ). with an inferred rotation speed 
that falls rapidly outside 300 arcseconds, accompanied by a 
rise in velocity dispersion. As previously noted, such kine- 
matics were not predicted by any of the simpler scenarios 
for SO formation, therefore motivating this further investi- 
gation. 

One immediate clue is offered by the only significant 
difference between the results of the two analyses presented 
in Figure (5] Specifically, while th e radial dispersion exceeds 
the tangential component in the INoordermeer et al.l (|2008l ) 
analysis, the opposite is the case in the maximum likeli- 
hood fit. In fact, the ordering of dispersions was fixed in the 
first analysis, as their ratio was set by the epicyclic approx- 
imation, which forc es Or > o"^ unless the rota tion curve is 
very rapidly rising (|Binnev fc Tremaine|[l987l ). In the cur- 
rent analysis, we have left both components as free parame- 
ters, and it is telling us that we then find > ar, which is 
not what the physically-motivated epicyclic approximation 
would produce, suggesting that there is some basic flaw in 
the model. 

As noted above, the principal missing element in this 
model is the spheroidal component, which is consistent with 
the enhanced value of a^,. Specifically, a spheroid will con- 
tribute a relatively small number of PNe with a large dis- 
persion but zero mean velocity. Close to the minor axis of 
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Figure 3. Derived circular velocity and the components of the velocity dispersion versus radius for the model galaxy in the disk and 
in the bulge. The left panels show the result when the spheroid's Sersic index is fixed at its true value of 3.5, while the right panels 
show the results when it is left as a free parameter. Upper panels arc for a larger catalogue of 437 kinematic tracers; lower panels for a 
smaller sample of 136. The filled symbols are from the maximum-likelihood analysis, with vertical error bars indicating uncertainty and 
horizontal error bars showing the extent of each radial bin (with the point plotted at the median location for a kinematic tracer in that 
bin). 



the galaxy, the disk PNe will also have zero mean line-of- 
sight velocity, whereas on the major axis their distribution 
will be offset in velocity due to rotation. Combining the two 
distributions with different mean velocities will result in a 
larger incorrectly-inferred dispersion than combining them 
where their mean velocities are the same. Since the tangen- 
tial component projects mostly into the line-of-sight close to 
the major axis, its derived value would be more enhanced 
by this contamination than that of the radial component, 
consistent with the results in Figure [S] 

A very simple test of whether such contamination could 
be responsible for the unphysical results can be made by 
increasing the severity of the likelihood clipping, to try to 



remove the contaminants from the fit. Figure [6] shows the 
success of such a process, in which a more aggressive thresh- 
old resulted in 34 PNe being rejected. The components of the 
velocity dispersions are now ordered in the manner physi- 
cally expected for a disk population. Further, the bizarre 
behavior of the kinematics has entirely vanished, with the 
rotation velocity now remaining approximately constant out 
to large radii, and the dispersion remaining low, just as one 
would expect for a normal disk population. 

As further evidence that the cause of the contamination 
is the spheroid. Figure |4] highlights the locations of the PNe 
rejected in this iterative clipping. The PNe appear spread 
throughout the galaxy, and not flattened into the aspect ra- 
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Figure 6. As for Figure[5l but with the hkehhood chpping prob- 
abiUty threshold increased to 1.65cr. 

tio of the disk, indicating that they are probably not the 
result of a poor disk model. There is, however, some indica- 
tion of an asymmetry in rejected PNe between the two sides 
of the galaxy, which we return to in Section 14.41 

Clearly we need to model a spheroidal component as 
well as the disk, but it is also heartening to see that the 
likelihood rejection method does a respectable job of deal- 
ing with even such a high level of contamination, in which 
almost 20% of the PNe do not seem to be drawn from the 
assumed model. This result provides some confidence in the 
robustness of the adopted procedure. 

4.2 Spheroid — disk decomposition 

The first step toward modeling the spheroidal compo- 
nent is to decompose continuum images into disk and 
spheroidal compo nents, which we have done using GALFIT 
l|Peng et al.l 12002 ) to fit an exponential disk and a de Vau- 
couleurs law spheroid to deep images of NGC 1023. Fig- 
ure [7] shows the result of this process carried out on both B 
and R band images. The companion galaxy NGC 1023A is 
mildly apparent even in the raw images, but shows up very 
clearly in the residuals when the model is subtracted. The 
residual image also shows e vidence for the bar at t he cen- 
tre of this galaxy noted bv iDebattista et al.l l|2002l ): since 
this faint feature is relatively localized at small radii, and 
we are primarily interested in the balance between disk and 
spheroid hght at large radii, we do not attempt to model it 
any further. We could also have considered other models for 
the spheroid, such as incorporating distinct bulge and halo 
components, or picking a more general Sersic profile for the 
light distribution. However, as we have seen in Section[3l the 
results are relatively insensitive to such subtleties. Since the 
primary goal is just to model to a reasonable approximation 
the fraction of the light in the different components at differ- 
ent positions, there is little benefit in trying to distinguish 
between what are likely to be near-degenerate alternative 
fits. 

From the models in the two bands, we can calculate col- 
ors for the individual components, and in this case we find 
B — R — 1.62 for the disk and 1.61 for the spheroid. Thus 
there appears to be essentially no color difference between 
the components. We also find no significant difference be- 
tween the scale-lengths of the components in the two bands: 
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Figure 8. Map showing the probability, /, that a PN at any given 
location is drawn from the spheroidal component of NGC 1023. 
The positions of the detected PNe are superimposed. Bold squares 
identify NGC1023A PNe. The / = 0.5 contour is shown in black. 



the best-fit model disk has a scale length of 67 arcseconds 
for the B-band data and 66 arcseconds for the R-band; sim- 
ilarly, the spheroid has an effective radius of 24 arcseconds 
in the B-band and 21 arcseconds in the R-band. There is 
thus no evidence for color gradients within each component. 
This lack of color variation between and within components 
makes the translation from stellar continuum properties to 
the probability that a PN belongs to the spheroid or disk 
particularly simple in this the probability is just 

the ratio of the component to total fight at each point. The 
bulge to total light ratio is 0.31 in the R-band and 0.35 in 
the B-band. These values compare well with estimates from 
the cruder on e-dimensional b ulge-disk fight decomposition 
performed in (|Noordermeer et al. 2008), where the bulge-to- 
total fight ratio was found to be 0.36. The two-dimensional 
map of the fraction of fight in the spheroidal component, 
/, as a function of position on the sky is presented in Fig- 
ure [51 As discussed in Section 12.31 the value of / clearly 
varies strongly with azimuth, with the spheroid being the 
dominant component at all radii close to the minor axis, 
underlining the necessity of this full two-dimensional decom- 
position. 

4.3 The disk + spheroid model 

Having calculated the division between spheroid and disk 
light, we can now carry out the likelihood analysis for 
NGC 1023 incorporating this extra component, as set out 
in Section 12.21 the resulting best-fit kinematic parameters 
as a function of radius are shown in Figure [9] We now start 
to see the characteristic properties of a normal disk galaxy. 
Rotation dominates random motions in the disk at all radii, 
and the mean rotation speed stays approximately constant 
out to the last points shown. In this plot, we have also filled 
in the missing stellar rotation velocity at small radii from 
conventional absorption- line data, and it is clear that the 
two techniques agree extremely well where they overlap. The 
spheroid dispersion profile is also very well-behaved, showing 
little if any variation with radius (justifying the compromise 
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Figure 7. GALFIT analysis of NGC 1023. The top row is for a R-band image; the bottom row is for an B-band image. The data are 
obtained with the Isaac Newton Telescope on 19 95 December 25, using the prime focus camera. The 300 second exposures cover a field 
of view of 10.2' X 10.2' jNoordermeer et al.ll200^ ). The first column shows the original data, the second column presents the model, and 
the third column shows the result of subtracting the model from the data. The residual image has been normalized by the original data, 
so that the values range from for a perfect fit, to 1 or —1, where the data has not been fitted at all. 



in radial binning discussed in Section 12. 2p . It is also no- 
table that the approximately-constant spheroid dispersion 
is consistent with a value of which is what one would 

obtain for the simplest possible model of a singular isother- 
mal sphere potential, in which circular speed and velocity 
dispersion are related in this way. This consistency is reas- 
suring, as the fitting procedure in no way imposed it on the 
results. 

As one further enhancement, one might also expect 
some rotational velocity, Vs, in the spheroidal component, 
which might affect the results. Indeed, there is known to 
be a strong correlation between ellipticity and rotational 
speed in low-lumi nosity spheroidal systems like this bulge 
(see Figure 4.14 in lBinnev fc Tremaine|[l987l '). The GALFIT 
modeling shows that the spheroid in NGC 1023 has an el- 
lipticity of ~ 0.25, which translates into a predicted value 
of Vs/cTj, ~ 0.5. We have experimented with including rota- 
tion at this level in the spheroidal component by modifying 
the first term in Equation [T] but find that at this level the 
inclusion of rotation makes no significant difference to the 
remaining kinematic parameters. 

The only slight disappointment in the fitting process is 
that the ratio between the components of disk dispersion is 
somewhat less well defined than was the case in the disk-only 
model fit, which presumably reflects the impact of the extra 
spheroidal dispersion free parameter in this model. However, 
now that we clearly have a well-behaved normal disk system 



in this galaxy, it seems reasonable to reduce the number of 
free parameters by invoking the epicyclic approximation. In 
particular, for such a cold disk system with a flat rotation 
curve, we expect a^/ar = l/V^. With this additional con- 
straint, we obtain the kinematic parameters plotted in Fig- 
ure [10] The error bars on all parameters are duly reduced, 
and we now find that the kinematics of NGC 1023 are ex- 
actly what one would expect for a very normal disk galaxy, 
with V/<Jti> ~ 4.1 throughout the disk, similar to what one 
finds in the stellar component of a spiral galaxy. Thus, we 
seem to have found the explanation for the strange kinemat- 
ics originally inferred from these data for NGC 1023, and it 
is now revealed to be most likely a spiral galaxy that has 
simply been stripped of its gas. 



4.4 PNe objects rejections: the footprint of an 
ongoing merger 

However, the story does not quite finish there. Even with 
the full disk-|-spheroid kinematic model, 17 PNe are still 
likelihood-clipped at a threshold probability of 2.1a. In Fig- 
ure 111! we show the locations of the PNe in both velocity 
and right ascension, with the rejected objects highlighted. 
Since NGC 1023 lies at a position angle very close to 90 de- 
grees, the spatial coordinate is essentially the distance along 
the major axis, so this is a conventional position-velocity di- 
agram, with the usual antisymmetric signature of a rotating 
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Figure 9. Derived mean rotation speed, the components of the 
disk velocity dispersion, and the spheroid dispersion versus ra- 
dius for NGC 1023. The filled symbols are from the maximum- 
likelihood analysis, with vertical error bars indicating uncertainty 
and horizontal error bars showing the extent of each radial bin 
(with the point plotted at the median location for a PN in that 
bin. The open symbol s show rotation veloci ties derived from 
absorption-line data bv IPebattista et al 
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Figure 10. As for Figure[9l but with the ratio of disk dispersions 
components constrained by the epicyclic approximation. 



disk evident in the PNe. The rejected PNe, however, do not 
show such antisymmetry, with the vast majority located on 
one side of the galaxy. Such an asymmetric arrangement is 
clearly not consistent with errors arising from a poor choice 
for any axisymmetric element in the model, or from mis- 
identified unrelated background objects. 

A clue to their origin comes from considering the loca- 
tion of NGC 1023A in this plot, as the rejected PNe seem 
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Figure 11. Plot of radial velocity versus right ascension for the 
PNe in NGC 1023. Filled symbols show the objects rejected by 
the likelihood clipping in the full disk-|-spheroid model. Asterisks 
show the PNe attributed to NGC 1023A. 



almost all to form a continuous stream that passes through 
this companion galax y, as one might expect i f the systems 
are tidally interacting ICapaccioli et al.l l|l986l ). It therefore 
seems likely that these PNe lie in the tidal debris from this 
companion as it is stripped in an on-going minor merger. 



5 CONCLUSIONS 

We have presented a new method for analyzing the kinemat- 
ics of disk galaxies as derived from individual stellar tracers 
such as PNe. This hybrid technique bins data radially in the 
galaxy to maintain the maximum flexibility in the inferred 
kinematics, but uses a likelihood analysis within each bin 
so as to derive the maximum amount of information from 
the discrete data points. In addition, we use photometric 
decomposition of continuum images to assign a probability 
to each kinematic tracer of belonging to either the disk or 
the spheroidal component of the galaxy being modeled. 

Application of the method to simulated data shows that 
it reproduces most of the intrinsic dynamics of a galaxy even 
when the number of discrete kinematic data points is rela- 
tively modest. Application of this technique to NGC 1023 
has offered an explanation for the strange kinematics pre- 
viously inferred for this system. These peculiar properties 
can be entirely attributed to uncorrected contamination by 
the galaxy's spheroidal component, and when this element 
is properly modeled, the galaxy is revealed to have stellar 
kinematics entirely consistent with a normal spiral system, 
with a cold rotationally-supported disk and a hot spheroid 
with the expected velocity dispersion. This result favours a 
model in which NGC 1023 formed from a spiral via a sim- 
ple gas stripping process or secular evolution, rather than 
through a more disruptive minor merger. 

One benefit of the likelihood analysis made clear by this 
application is that it is possible to "likelihood clip" , to iden- 
tify objects that do not seem to be drawn from the model, 
in a reasonably robust and objective manner. Comparison 
of Figures [TO] and [5] shows that such clipping can deal quite 
well with even the rather extreme case where only the disk is 
included in the model. Where a more realistic disk-|-spheroid 
model is adopted, it seems to do a good job of identifying 
features like stellar streams, adding to our understanding of 
the dynamical properties of such systems. 
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Clearly, though, although the detailed dynamics of this 
galaxy arc interesting, one cannot infer any general conclu- 
sions about the formation of SO galaxies from a single object. 
Therefore, the next step in this project will involve applying 
this analysis technique to observations of a larger sample of 
SO galaxies whose PNe kinematics have been observed with 
PN.S. By obtaining a measure of the stellar kinematics of 
SOs in regions spanning a range of galaxy densities, we will 
be able to see if they all have the stripped-spiral properties 
of NGC 1023, and hence whether there is a single route to 
SO formation irrespective of environment. 
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